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ABSTRACT 

The bispectrum B(ki, k%, k^), the three-point function of density fluctuations in Fourier space, 
is the lowest order statistic that carries information about the spatial coherence of large-scale 
structures. For Gaussian initial conditions, when the density fluctuation amplitude is small 
(8 <C 1), tree-level (leading order) perturbation theory predicts a characteristic dependence of 
the bispectrum on the shape of the triangle formed by the three wave vectors. This configuration 
dependence provides a signature of gravitational instability, and departures from it in galaxy 
catalogs can be interpreted as due to bias, that is, non-gravitational effects. On the other hand, 
Af-body simulations indicate that the reduced three-point function becomes relatively shape- 
independent in the strongly non- linear regime (5 > 1). 

In order to understand this non-linear transition and assess the domain of reliability of shape- 
dependence as a probe of bias, we calculate the one-loop (next-to- leading order) corrections 
to the bispectrum in perturbation theory. We compare these results with measurements in 
numerical simulations with scale-free and Cold Dark Matter initial power spectra. We find 
that the one-loop corrections account very well for the departures from the tree-level results 
measured in numerical simulations on weakly non-linear scales (5 < 1). In this regime, the reduced 
bispectrum qualitatively retains its tree-level shape, but the amplitude can change significantly. 
At smaller scales (6 > 1), the reduced bispectrum in the simulations starts to flatten, an effect 
which can be partially understood from the one-loop results. In the strong clustering regime, 
where perturbation theory breaks down entirely, the simulation results confirm that the reduced 
bispectrum has almost no dependence on triangle shape, in rough agreement with the hierarchical 
ansatz. 



Subject headings: large-scale structure of universe; methods: numerical; methods: statistical 
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1. Introduction 

The growth of cosmological density fluctuations in perturbation theory (PT) is becoming a mature, well- 
understood subject, with techniques established that in principle allow calculations to arbitrary order (e.g., 
Goroff et al. 1986; Jain & Bertschinger 1994). On large scales, where the rms density fluctuations are small, 
tree-level (leading order) PT gives the first non- vanishing contribution to statistical averages, and it has been 
used to understand the generation of higher order correlations in gravitational instability (e.g., Peebles 1980; 
Fry 1984; Juszkiewicz, Bouchet & Colombi 1993; Bernardeau 1992, 1994). Comparison with fully non-linear 
numerical simulations has shown this perturbative approach to be very successful (Juszkiewicz et al. 1993, 
1995; Lucchin et al. 1994; Bernardeau 1994; Fry 1994a; Lokas et al. 1995; Gaztahaga & Baugh 1995; Baugh, 
Gaztahaga & Efstathiou 1995). 

On smaller scales, where the density fluctuation amplitude approaches or exceeds unity, loop (next-to- 
leading and higher order) corrections to the tree-level PT results should become important. The question 
then arises of whether our understanding of clustering can be extended from large scales further into the 
non-linear regime by using one-loop PT. In previous papers, we considered one-loop corrections to one-point 
cumulants of unsmoothed fields (Scoccimarro & Frieman 1996), the power spectrum, variance, and two-point 
correlation function (Scoccimarro & Frieman 1996b; see also Makino, Sasaki, & Suto 1992; Lokas et al. 1996), 
and the bispectrum and skewness including smoothing effects (Scoccimarro 1997). 

In this paper, we present one-loop perturbative bispectra for a variety of initial power spectra, and we 
compare them with results of numerical simulations. The bispectrum B(k±, &2, k%), the Fourier transform 
of the connected three-point correlation function of the density field perturbations, is of interest for several 
reasons. For Gaussian initial conditions, the connected N— point correlation functions vanish in linear theory 
for N > 2. The bispectrum is therefore intrinsically non-linear, and it is the lowest-order statistic with this 
property. It is also the lowest-order statistic that carries information about the spatial coherence of the 
density and velocity fields; by contrast, the density power spectrum P(k) ~ (|<5(fc)| 2 ) is independent of 
phase correlations. 

In addition, the dependence of the tree-level bispectrum on the shape of the triangle formed by the 
three wave vectors k,\,k,2, is a characteristic signature of gravitational instability. The degree to which the 
observed galaxy distribution exhibits this predicted configuration dependence provides an independent probe 
of bias, that is, of the relation between the galaxy and mass density distributions (Fry 1994b). In the galaxy 
data that have been studied to date, the Shane- Wirt anen (Lick) angular catalog, the reduced bispectrum 
was not found to exhibit the expected tree-level dependence on configuration shape (Fry & Seldner 1982), 
and this can be interpreted as a possible indication that these galaxies are biased relative to the mass (Fry 
1994b). However, numerical simulations indicate that the reduced bispectrum becomes much less dependent 
on configuration (Fry, Melott, & Shandarin 1993, hereafter FMS; Fry, Melott, & Shandarin 1995) in the 
highly non-linear regime. Thus, in order to rigorously apply this test, i.e., to use the bispectrum shape as 
a quantitative probe of bias, we must assess the reliability of the tree-level predictions and quantify where 
and how they break down. One-loop PT, in conjunction with A^-body simulations, provides a framework for 
accomplishing this. 

Study of the one-loop bispectrum is also of practical value for the analysis of galaxy surveys. For 
scale- free initial power spectra, P(k) oc k n , measurement uncertainties in the reduced bispectrum scale as 
(k ne /k)( n+ V/ 2 (FMS), where k is an inverse wavelength and k n i is the scale of nonlinearity (see eq. |27|] 
below). Therefore, measuring the reduced bispectrum is difficult on large scales (k <C k n e). For better 
accuracy, we would like to measure B on scales with k as large as possible, at least not much smaller than 
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k n f. Since the tree-level PT result is not expected to be accurate in the strong clustering regime, k > k n e, we 
would like to know the size and nature of the corrections to the tree-level prediction on scales k comparable 
to k n £. This is precisely the domain where one-loop PT can be very helpful, and can be additionally checked 
against numerical simulations which arc reliable at these scales. 

The paper is organized as follows. In Section 2, we briefly review PT solutions and the loop expansion 
of the power spectrum and bispectrum. In Section 3, we apply PT to power-law initial power spectra and 
give some analytic results. Section 4 describes the numerical simulations analyzed in this work and Section 5 
compares them to the one-loop perturbative results on the power spectrum and bispectrum. Section 6 
contains a final discussion. Finally, Appendix A describes the method used to measure the bispectrum in 
numerical simulations, and Appendix B reviews the general framework of PT, including the extension to 
arbitrary cosmological parameters £1 and A. 



2. Perturbation Theory 

We work in transform space with the Fourier amplitude of the density contrast S(x) = [p(x, i) — p]/p, 
defined such that 

To linear order, S(k,t) = a(i)8i(k), where a(t) is the cosmic scale factor (normalized to a(t ) = 1 today), 
and 6i(k) denotes the linear density fluctuation amplitude (at the present epoch with the scale factor 
normalization above). We consider an Einstein-de Sitter Universe, with density parameter f2 = 1, in which 
case the PT expansion for the density contrast can be written 

oo 

5(fc,i) = <Ufc)- (2) 

n=l 

Modeling the matter as pressureless non-relativistic 'dust', an appropriate description for cold dark matter 
before shell crossing, the fluid equations of motion (see Appendix B.l) determine S n (k) in terms of the linear 
fluctuations, 

S n (k)= J <!%... J d 3 q n S D (k- qi q n )F^(q 1 ,...,q n )5 1 (q 1 )---6 1 (q n ), (3) 

where the are dimensionless, symmetric, scalar functions of the wave vectors {q l , . . . , q n } (Goroff et al. 
1986; Jain & Bertschinger 1994) and 6u is the Dirac ^-function. The recursion relations from which these 
kernels can be derived are provided for reference in Appendix B. 

A systematic framework for calculating correlations of cosmological fields in PT has been formulated 
using diagrammatic techniques (Fry 1984; Goroff et al. 1986; Wise 1988; Scoccimarro & Frieman 1996). 
From this point of view, leading order PT for the statistical quantities of interest corresponds to tree graphs, 
next-to-leading order PT contributions can be described in terms of one-loop graphs, etc. 

The simplest statistic of interest is the power spectrum P(k), the second moment of the Fourier ampli- 
tude of the density contrast, defined by 

(6(k)6(k')) =Su(k + k')P(k). (4) 

By statistical isotropy, the power spectrum depends only on the magnitude of k. To tree level (linear theory), 
the power spectrum keeps its shape and simply grows by an overall amplitude. One-loop (first non-linear) 
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corrections introduce coupling between different Fourier modes and increase or decrease the growth rates 
relative to linear theory, depending on the amount of small-scale power in the initial spectrum (Klypin & 
Melott 1992; Makino et al. 1992; Lokas et al. 1996; Scoccimarro & Frieman 1996b). We can write the loop 
expansion for P(k) up to one-loop corrections as (henceforth we suppress the implicit time dependence) 

P{k) = P {0 \k) + P {1 \k) + ■ ■ ■ . (5) 

The superscript (n) denotes an n-loop contribution. The tree-level (0-loop) contribution is just the linear 
spectrum, with 

a 2 (S^S^k 1 )) =6 B (k + k')pM(k), (6) 
and the one-loop contribution consists of two terms, 

P (1) (fc)-P 22 (fc)+P 13 (fc), (7) 

where 

P 22 (k) = 2|[F 2 (s) (fc-q,q)] 2 p(°)(|fe-q|)p( )( (Z )A, (8) 

Pl3(k) = 6jFt } (k, q ,-q)P^(k)P^( q )d 3 q . (9) 

Here Pij denotes the amplitude corresponding to the contribution (Si(k)Sj(k)) to the power spectrum. We 
have assumed Gaussian initial conditions, for which Pij vanishes if i + j is odd. 

The third moment in the Fourier domain gives the bispectrum, B(k\,k 2 , k 3 ), defined by 

(5(k 1 )5(k 2 )5(k 3 )) = S D (k 1 + k 2 + k 3 )B(k 1 ,k 2 ,k 3 ). (10) 

The Dirac S- function ensures that the bispectrum is defined only for configurations that form closed triangles, 
^2 ki — 0. For Gaussian initial conditions, the first non-vanishing contribution to the connected n-point 
correlation function requires PT to order n — 1 (Fry 1984). The loop expansion for the bispectrum reads 
(Scoccimarro 1997): 

B(k u k 2 , h) = B<°) (fci, hi, k 3 ) + (fci, k 2 , fc 3 ) + • • • • (11) 

The tree-level term is 

B (0) (*i, k 2 , k 3 ) = 2 P^(h) P^{k 2 ) F 2 (s) (fcx, k 2 ) + 2 P(°)(fe2) P(°)(fc 3 ) P 2 (s) (fc 2 , fe 3 ) 

+2p( )(fc 3 )p(°)(fc 1 )P 2 (s) (fe 3 ,fc 1 ), (12) 

where the second-order kernel for the Einstein-de Sitter model is 

2F«(fe i) kj) = y + k ■ fc,(| + I) + \{k t ■ kjf ; (13) 

hats denote unit vectors (Fry 1984). For Q ^ 1, and vanishing cosmological constant, the factors 10/7 and 
4/7 become 1 + k and 1 — k, where k ~ |fi~ 2 / 63 depends only very weakly on f2 (Bouchet et al. 1995, 
Hivon et al. 1995). A discussion of the dependence of PT kernels on cosmological parameters is presented 
in Appendix B.3. 

The one-loop bispectrum comprises several terms, 

B^\h, k 2 ,k 3 ) = B 222 (h, k 2 , k 3 ) + S| 21 (fci, fc 2 , fc 3 ) + B^ih, k 2 ,k 3 ) + B 4 u(h,k 2 ,k 3 ), (14) 
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with: 



B 222 = 8 J d 3 qP^(q)F^(-q,q+k 1 )P^(\q + k 1 \)F^(-q-k 1 ,q-k 2 ) 

xP^(\q-k 2 \)F^(k 2 -q,q), (15) 
BL = 6P {a Hk3)Jd 3 qP^( q )Ft\-q,q-k 2 ,-k 3 )P^(\q-k 2 \) 

xF% (q,k 2 — q) + 5 permutations, (16) 
BUi = 6 P(°) (k 2 ) P(°) (* 3 ) F 2 (s) (fe 2 , fe 3 ) | d\ (g) F 3 (s) (fc 3 , q, -q) 

+ 5 permutations, (17) 
P411 = 12 P (0) (fc 2 ) P (0) (A*) / d 3 qP^(q)Fi s) (q,-q,-k 2 ,-k 3 ) 

+ 2 permutations. (18) 

Again, the subscript ijk denotes a contribution of order (SiSjSk) to the bispectrum. The reader is referred 
to Scoccimarro (1997), Fig. 3, for a diagrammatic representation of these terms. 

The reduced bispectrum, or hierarchical three-point amplitude Q is defined as 
n(h 1 , n_ B(k u k 2 ,k 3 ) 

W U 2 ' 3j ~ P(fci) P(fca) + P(fc 2 ) P(fc 3 ) + P(fc 3 ) P(fci) ' UJj 



The loop expansion of the numerator and denominator yields: 

^_ BW(k u k 2 ,k 3 ) + B^(k 1 , k 2 ,k 3 ) 



x(°Hk 1 ,k 2 ,k 3 ) + zW(k 1 ,k 2 ,k 3 ) + 



(20) 



where 



E<°> (*!,**,**) = pW(fci)P (0) (fe) + p(°)(A :2 )p( )(fc 3 ) + p( )(A : 3)P (0) (fci), (21) 
E (1) (fci,/fc 2 ,/c 3 ) = P (0) (fci)P (1) (/c 2 ) + 5 permutations. (22) 



The loop expansion of Q = + Q^ 1 ' + ■ ■ ■ gives 



)(0) = B^(h,k 2 ,k 3 ) 

^){k u k 2 ,k 3 y 



(23) 



w - E (o) S (o) • ^ 

Note that depends on the normalization of the linear power spectrum, and its amplitude increases with 
time evolution. On the other hand, from equations (|l2|), (pl|), and ( p3| ) it follows that is independent 
of time and normalization (Fry 1984). Furthermore, for scale-free initial conditions, P(°\k) <x fc™, 
is also independent of overall scale. For the particular case of equilateral configurations (fci — k 2 = k 3 
and ki ■ kj = —0.5 for all pairs), is independent of spectral index as well, Q^q — 4/7. In general, 
for scale- free initial power spectra, depends on configuration shape through, e.g., the ratio k\/k 2 and 
the angle 9 defined by k\ ■ k 2 = cos6. This configuration dependence of reflects the anisotropy of 
structures and flows generated by gravitational instability. From equations (p3|), (|2l|), and (021) , it follows 
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that would be independent of configuration shape if F% (ki,kj) (see eq. 13 ) were a constant. The 
configuration dependence of F% (ki, kj) implied by equation ( |l3| ) has two sources: the term linear in ki ■ kj 
comes from the gradients of the density field in the direction of the flow, whereas the term quadratic in ki ■ kj 
represents the gradients of the velocity field in the direction of the flow (Scoccimarro 1997). Thus, is 
enhanced if the wave vectors are collinear (6 — 0, 7r), which reflects the fact that large-scale flows generated 
by gravitational instability are mostly parallel to density gradients (see, e.g., dotted lines in Fig. 1). This 
physical interpretation provides some insight into what is expected to happen as the transition to the non- 
linear regime is made. As long as the evolution of structures is dominated by large-scale motions, the shape 
dependence of Q should remain qualitatively the same as at tree level. In fact, for spectra dominated by large 
scale power, the anisotropy of structures is amplified by the "pancaking" process described by the Zel'dovich 
(1970) approximation (Coles et al. 1993; Melott & Shandarin 1993) and this would lead to an enhancement 
of the configuration dependence of Q. On the other hand, when substantial velocity dispersion develops 
on small scales due to virialization, the interpretation above suggests that one should see a flattening in Q: 
non-collinear configurations become more probable, due to the loss of coherence of structures and flows (i.e., 
the gradient terms in eq. |l3| are smoothed out due to random motions.) 

To characterize the degree of non-linear evolution when including one-loop corrections to the power 
spectrum and bispectrum, it is convenient to define a physical scale from the linear power spectrum. One 
such scale is the correlation length Rq, the scale on which the smoothed linear variance is unity, af(Ro) = 1. 
The variance is defined by 

a 2 e (R)= f d 3 kP {0 \k)W 2 (kR), (25) 



where W(x) is the Fourier transform of a window function (usually a top-hat or Gaussian) of characteristic 
scale R. For scale-free initial power spectra, P( \k) — Aa 2 k n , the variance scales as c 2 (R) — (R/ Rq)' 1 -- 71 ^. 
For Gaussian smoothing, W(kR) — exp(— k 2 R 2 /2), the linear correlation length satisfies 

R^^2nAa 2 r( 7 -^). (26) 



In the 128 3 PM simulations described below the epoch of evolution is labeled by k n £, the wave number that 
is on the threshold of going nonlinear as determined in linear theory, defined by 



4tt 



d ^P^(k) = J -^Aa 2 k:l^l, (27) 



where kf = 2ir/L is the fundamental mode of the simulation box of side L. Equations (|26|) and ( |27| ) imply 
that k n( R Q ^ r[(n + 5)/2]. 



3. One-Loop Results for the Bispectrum 

For scale-free initial power spectra, the one-loop integrals in equations (|l5|)-([l^) can be calculated 
analytically in the range —3 < n < —1 by using dimensional regularization (Scoccimarro 1997). In this 
spectral range, the resulting bispectrum obeys self-similarity and, based on previous results for the power 
spectrum (Scoccimarro & Frieman 1996b), the one-loop calculations are expected to give a good description 
of the transition to the nonlinear regime. 

Due to statistical homogeneity and isotropy, the bispectrum B(k±, ki, ^3) depends on time, the mag- 
nitudes k\, &2, and the angle 9 (ki ■ ki = cos#). In order to display the analytic results, however, 
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it is more convenient to trade the variable 8 for the third side of the triangle, fc 3 = \ki + fca|- Let 
SW(fci,fe2,fe 3 ) = A 3 a 6 7r 3 &W(fci,fc 2 ,fc 3 ), with fci + fc 2 + fc 3 = 0. Then, for n = -2 it follows (Scocci- 
marro 1997): 



30279 
34496 fci 3 

233 fci 6 
8624 fc 2 4 fc 3 5 

23 fci 5 
103488 fc 2 4 fc 3 4 

5311 fci 2 
34496 fc 2 2 fc 3 3 "* 



2635 ki 



37313 fci 



38431 



51744 k 2 b 206976 fc 2 



68992 fci fc 2 



16517fci 5 

* — F 4 

362208 fc 2 3 fc 3 5 
9791 fci 4 
206976 fc 2 3 fc 3 4 
42983 fci 
362208 fc 2 fc 3 3 "* 



7392 fc 2 2 fc 3 5 



197 fci 78691 fci 3 

275968 fc 2 fc 3 5 
703 fci 3 , 19867 fci 2 
206976 fc 2 k 3 4 
28393 
19712 fci fc 2 fc 3 



68992 fc 2 2 fc 3 4 
131 fci 



3696 fc 2 2 fc 3 2 



53973 fci 7 108685 fci fc 2 59599 fci 3 



1931776 fc 2 5 fc 3 5 
permutations. 



181104 fc 3 5 



362208 fc 2 3 fc 3 3 



(28) 



A simple result can be obtained for equilateral configurations. Given that the one-loop power spectrum 
for n = —2 can be written as pW(fc) = A 2 a 4 557r 3 /(98fc) (Scoccimarro & Frieman 1996b), the hierarchical 
amplitude for equilateral configurations at the one-loop level is: 



Q EQ (n=-2) = - 



1426697 



tt 3/2 kR = 0.57 + 2.06 fci? - 



For n 



3863552 

-1.5, the corresponding result reads: 

Q EQ (n = -1.5) = 0.57 + 1.32 (fci? 



,3/2 



(29) 



(30) 



For other spectral indices in the range —3 < n < — 1, the one- loop bispectrum can be expressed in terms 
of hypergeometric functions (Scoccimarro 1997). On the other hand, when n > — 1, one-loop PT leads to 
ultraviolet (fc — > oo) divergences that must be regulated by the introduction of a cutoff scale. In this case we 
take the initial power spectrum to be (fc) = A a 2 fc™ for e < fc < fc c and zero otherwise. For convenience 
in comparison with the n = — 1, 0, 1 numerical simulations described below, in these cases we take fci = 1, 
fc 2 = 1/2, e = 1/16, and fc c = 4. The integration of equations (|l^)-(0) is then done numerically. Given 
the complexity of the calculations involved, it is desirable to verify that the numerical integration code is 
correct. For this reason, we have written two completely independent codes of numerical integration, one 
based on Romberg integration, the other using Gaussian adaptive integration. The results of both codes 
agree with each other very well over the whole range of spectral indices and cutoff parameters considered. 
They also agree with the analytic result of equation ( p^ ) for the case n = —2, in the limit that the spectral 
cutoffs are removed. The one-loop bispectrum for cold dark matter (CDM) initial spectra does not present 
any additional complications for the numerical evaluation and is calculated using the same program. Since 
in this case w fc~ 3 at large fc, there are no ultraviolet divergences for this spectrum. Typically, the 
numerical evaluation of the one-loop bispectrum requires a few hours in one processor of a Silicon Graphics 
Power Challenge or DEC- Alpha 2100 workstation. 
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4. Numerical Simulations 

4.1. Scale- Free Simulations 

We compare our perturbative calculations to numerical results from two sets of scale-free simulations 
with Gaussian initial conditions. 

The first is an ensemble of simulations with initial spectral indices n = — 2, — 1, and +1, performed by 
Melott & Shandarin (1993) with the Particle-Mesh (PM) code of Melott (1986). These simulations involve 
^Vpar = 128 3 particles and a 128 3 staggered mesh. The force on a cell is a result of differencing the potential 
in its 8 neighbors. These PM simulations have about twice the usual resolution by using a staggered mesh 
scheme (Melott 1986; Melott, Weinberg & Gott 1988). Models with n = +1 and —1 have an initial amplitude 
such that the rms fluctuation averaged over one mesh cell is (Ap/p) — 0.05. Models with n = and —2 
have the slightly larger initial amplitude (Ap/p)o = 0.25. The power spectrum and bispectrum data we 
use here were measured in these simulations respectively by Melott & Shandarin (1993) and by FMS. Four 
independent realizations were generated for each value of n, using the same sets of random number seeds. 
Besides improving the signal, averaging over four realizations allows us to estimate fairly the uncertainties in 
our results. For these simulations, the epoch of evolution is labeled by k n i, the wave number that is on the 
threshold of going nonlinear as determined in linear theory given in equation ( p7f ) . The simulation output 
times are characterized by k n £ = 64, 32, 16, 8, and 4, or k ny /k n e = 1, 2, 4, 8, 16, where k ny = 64 is the 
Nyquist frequency of the simulation (wave numbers given in units of the fundamental mode of the simulation 
box of side L, kf = 2ir/L). 

We have also analyzed two simulations, one with n = —2 (already used in Lokas et al. 1996) and the 
other one with n = —1.5, done with a vectorized PM code (Moutarde at al. 1991) modified to run in parallel 
on several processors of a CRAY-98 (Hivon 1995). They involve 256 3 particles and use a 256 3 (unstaggered) 
mesh to compute the forces. The initial conditions were set by using the Zel'dovich approximation on a 
"glass" (see, e.g., White 1994), and the initial power spectrum is given by P(k) = A 2 (n) (k / k ny ) n / 256 3 , 
where k ny is the Nyquist frequency of the particles in units of the fundamental mode, and A{— 2) = 0.2, 
A(— 1.5) = 1/V128. For n = —2, we have analyzed outputs with a — 8, 11.31, 16, 22.63, whereas for n = —1.5 
we have analyzed outputs with a — 22.63,32,45.25,64, where both simulations start at a = 1. To avoid 
spurious effects (Melott et al. 1988; Kauffman & Melott 1992; Colombi, Bouchet & Schaeffer 1994, 1995), we 
only consider scales k that satisfy the requirement 4 < k\/kf < k ny /2, where kf = 2tt/L is the fundamental 
mode of the simulation box. 



4.2. CDM Simulations 

The CDM simulation we analyzed was done by Couchman, Thomas & Pearce (1995) with an adaptive 
P 3 M (Particle-Particle-Particle-Mesh) code and involves 128 3 particles in a box of length 100 hr 1 Mpc (h = 
ifo/lOOkms -1 Mpc -1 , where Ho is the Hubble constant). These simulation data are publicly available 
through the Hydra Consortium Web page ( ^ittp: / /coho. astro. uwo.ca/pub/consort. html ). They correspond to 



an n = 1 model, with linear CDM power spectrum characterized by a shape parameter T = Qh = 0.25, in 
approximate agreement with the observed galaxy power spectrum on large scales (e.g., Peacock & Dodds 
1994). The initial conditions were set by using the Zel'dovich approximation on a grid. We have analyzed 
output times at which the linear rms density fluctuation amplitude in top-hat spheres of radius R = 8/t" 1 Mpc 
is given by <r$ — 0.2057, 0.3291, 0.64. These times correspond respectively to scale factors a = 0.3214, 0.5143, 
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1, where a = 0.02 initially. For the measurements, we have only considered scales in the range 4 < k/kf < 50. 
In particular, the k\/k% = 2 configurations shown in the hgures below correspond to k±/kf — 15, 30, 40. 

5. One-Loop Perturbation Theory vs. Numerical Simulations 

We now compare the one-loop perturbative predictions with the TV-body results for scale-free and CDM 
power spectra. Figure [j] presents our main results for the n = — 2 spectrum. On the top left panel, we show 
the contribution per logarithmic wave number interval to the variance, A(fc) = 4irk 3 P(k), as a function of 
scale. In linear theory, from equation (p7j), A' ) (fc) = (n + 3)(fc/fc„^) n+3 . The symbols correspond to the 256 3 
A-body results averaged over the four different time outputs, assuming self-similar evolution. The error bars 
in this plot are calculated from the dispersion in this averaging procedure (see Appendix A. 2 for details). 
In the same plot, we include the linear theory extrapolation, A^°'(k;n = —2) = (k/k n i), the one-loop PT 
prediction, and the phenomenological A-body fitting formulae proposed by Jain, Mo & White (1995, JMW) 
and by Peacock & Dodds (1996, PD), based on earlier work by Hamilton et al. (1991). Note how well the 
one-loop analytic result, 

A(kR ) = (2/y/n) kR (l + j^^' 2 = 1-128 kR (l + 1.562 fci? ) , (31) 

describes the numerical simulation measurements and fitting formulae up to scales where A ss 30, where 
the power spectrum goes over to the stable clustering regime. This is remarkable, given the simplicity of 
equation ([II]) when compared to the JMW and PD fitting formulae. 

In the other three panels in Figure |l|, we show results for the hierarchical amplitude Q for triangle 
configurations with k\/k% — 2, as a function of the angle between k\ and for three different scales 
corresponding to A(fci) = 0.71,1.95,3.92. Filled squares denote averages over the four simulation output 
times as mentioned above for the 256 3 PM code. Filled triangles denote results taken from the 128 3 PM 
simulation. In this case, the displayed values correspond to the average over the four different realizations, 
using only a single output time for each one. The error bars correspond to the dispersion over the mea- 
surements (see Appendix A). They should be more reliable than those estimated from the 256 3 simulations, 
which may also reflect departures from self-similarity. 

We note the clear departure of the A-body results from the tree-level PT prediction equation (^3|) at 
k\Ro = 0.45 and the very good agreement when the one-loop correction (eq. [Q) is included. The net effect 
of non-linear evolution at this stage is to partially increase the configuration dependence of Q, as expected 
from the fact that for this initial spectrum the evolution on weakly non-linear scales is still dominated by 
large-scale coherent motions. When A(fci) > 1, it is not justified to expand the denominator in equation 
( p0| ) to get equation ( pi| ) . Since the one- loop power spectrum agrees with that in the numerical simulations 
down to scales where the one- loop correction dominates over the tree-level contribution, when A(fei) > 1 wc 
use the full expression in equation ( |20| ) , denoted as "one-loop (s) ," which saturates at large kR due to self- 
similarity (Scoccimarro 1997). The two bottom panels in Figure [l] illustrate this situation for A(fci) = 1.95, 
3.92, where the simulation results show a gentle flattening of Q(9) as we approach more non-linear scales. 

In Figure ^ we present a similar set of plots corresponding to the n = —1.5 spectrum. The one- loop 
power spectrum is given by (Scoccimarro & Frieman 1996b): 

A(fci? ) = 1-632 (kR a ) 3/2 (l + 0.391 (A:i? ) 3/2 ) , (32) 

which (see top left panel) again shows very good agreement with the numerical simulations and the JMW and 
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PD fitting formulae up to scales where A(fc) ~ 10. In the remaining three panels, we show the hierarchical 
amplitude Q for stages of non-linear evolution similar to those of the corresponding panels in Figure |l|. We 
see the same trend with scale, namely, a departure from the tree-level PT prediction and then a hint of 
decrease in the configuration dependence of Q. Note that, in this case, the one-loop corrections to Q are 
smaller than in the n = —2 case, both for the power spectrum and bispectrum. In Figure |^, we show Q 
approaching the strong clustering regime, A(k) 3> 1, for both n = —2 (top panels) and n = —1.5 (bottom 
panels). These simulation results confirm the flattening of Q seen in previous work at small scales (FMS). 
Interestingly, the 9 ~ configurations (which correspond to more non-linear scales) seem to flatten before 
the 9 ~ 7r configurations. A similar effect is apparent in Figures |l| and ||: the loop corrections to the tree- level 
result are consistently larger for 9 ~ configurations. 

We also see from Figures [l], |[ and || that the flattening of Q happens first at intermediate angles, and 
then spreads to smaller and larger values of 9. This effect can be interpreted as follows: as explained in § 2 
(see discussion after eq. ^JJ), on weakly nonlinear scales, before shell-crossing, large-scale flows are mostly 
parallel to density gradients, an effect which favors collinear configurations (9 — 0, n). On smaller, more 
non-linear scales, the previrialization associated with shell-crossing leads to a "randomization" of gradients, 
i.e., configurations which do not pick out a preferred direction are given relatively more weight. This helps 
to explain why the flattening of Q first develops at intermediate values of 9. 

Figures ^ and |H show the results for n = — 1, 0, 1 spectra from the 128 3 PM simulations. As the spectral 
index n increases, the relative uncertainties in Q increase (see FMS, or Appendix A. 2, eq. [ A 1 7] ] ) , and it 
becomes more difficult to measure Q accurately in the simulations. Also, one-loop PT does not work well for 
these spectra, since ultraviolet divergences must be regulated by introducing a small scale cutoff, k c , which 
violates self-similarity. If the one-loop correction to the power spectrum were plotted in the top left panel 
in Figure |], we would have different predictions for different values of k c Ro (see Scoccimarro & Frieman 
1996b). The solid lines in Fig. |||show that by choosing the small-scale PT cutoff as the Nyquist frequency, 
it is not possible to match the results of the numerical simulation measurements. On the other hand, we 
see that the non-linear corrections found in the iV-body data are very small for n = — 1, for both the power 
spectrum and Q, and tree- level PT does well even on scales for which Awl. 

A similar situation is shown in the top panels of Figure || for the n = case. For n = 1 (bottom panels 
of Figure |^), the simulation data are quite noisy, which makes it difficult to reach any conclusion. Note, 
however, that the measured Q tends to be systematically above the tree- level prediction for n = +1, and 
actually agrees quite well with the tree- level predictions for n = (solid lines). This might be explained 
by the following argument. The initial spectrum of an n = +1 simulation is given by P(k) cx k up to 
the Nyquist frequency of the particles. For k > k ny , however, the initial ./V-body spectrum is generally 
white noise, n ff = 0, or something close to it, depending on the initial setup of the particle positions (e.g., 
Juszkiewicz et al. 1993; Colombi, Bouchet & Hernquist 1996). Usually, the details of the shape of the initial 
power spectrum at these small scales are unimportant, since power is expected to cascade from large to small 
scales, eventually establishing nearly correct small-scale behavior regardless of the initial state at such scales 
(e.g., Beacom et al. 1991; Little et al. 1991; Melott & Shandarin 1993; Bagla & Padmanabhan 1997). Since 
ann = +1 simulation does not have much relative power at large scales, however, the relaxation time for the 
cascade may be long enough for the initial conditions to affect three-point statistics such as the bispectrum 
at relatively late output times. Moreover, this effect may not be manifest in the power spectrum, which 
is insensitive to phase correlations. In this respect, the analysis of Fry, Melott & Shandarin (1992) and 
FMS suggests that, once the system has relaxed far enough into the nonlinear regime, the bispectrum is not 
significantly affected by the small-scale behavior of the initial conditions near the Nyquist frequency, even 
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for n = +1. (Their conclusions are based mainly on the analysis of equilateral configurations, however, and 
it would be interesting to extend this study to other configurations.) 

Figure [| shows the results from the CDM simulation for erg — 0.2057, the earliest output we analyzed 
in this case. For the linear CDM spectrum, we use the BBKS transfer function (Bardeen et al. 1986). For 
this weakly nonlinear output, where the amplitude of the power spectrum is not large compared to the white 
noise level, we have not corrected for the discrete nature of particles. Since the initial conditions are set by 
Zel'dovich displacements from a grid, Poisson noise is not a good model at early times (see, e.g., Melott & 
Shandarin 1993; Baugh & Efstathiou 1994). In fact, the standard shot noise correction would make the non- 
linear power spectrum even smaller than the linear spectrum. The excellent agreement of the uncorrected 
A-body power spectrum with one-loop PT and the PD fitting formula seems to indicate that this is the 
correct procedure. For the hierarchical amplitude Q, the effect of a full correction for discreteness would be 
smaller than on P(k): the corresponding values of Q in Figure^ would be somewhat higher than the results 
shown, in better agreement with the one-loop calculation for 9 = and 9 = tt (we do not plot this result for 
reasons of clarity.) All the other measurements in this paper have been corrected for discreteness. Overall, 
Figure ^ shows very good agreement between one-loop PT and the numerical simulation measurements, even 
when the deviations from the tree-level PT predictions are dramatic, as can been seen on the bottom right 
panel, corresponding to A(fci) = 1.05. The increase in configuration dependence of Q with A(fci) is more 
important than for the n = —2 case, consistent with what one would expect from the discussion in § 2 and 
the effective spectral indices n c ff(fc) = dlnP' ^/(ilnfc of the scales considered (displayed in Fig. ||). Note that 
the error bars on the plots probably underestimate the true errors. Having access to only one realization, 
and without the possibility of using self-similarity (and thus different output times) as a test on the accuracy 
of the results because the CDM spectrum is not scale-free, we have estimated the error bars from the number 
of independent modes in fc-space contributing to a given configuration (see Appendix A for details.) 

Figure shows the set of plots corresponding to the next output time analyzed in the CDM simulation. 
We see that the logarithmic variance contribution A at one-loop agrees very well with the TV-body measure- 
ments up to A(fc) ss 6. For the hierarchical amplitude Q, given that we have only one realization, we could 
only make accurate measurements on scales for which A(fci) > 1, so we use equation ( p0| ) for the one-loop 
prediction. We see very good agreement between predictions and measurements for configurations close to 
collincar (9 = 0, tt), and a progressive flattening of Q as we move to smaller scales, in agreement with the 
results for n = — 2 and n = —1.5. The latest CDM output, shown in Figure ||, illustrates this behavior 
of Q further in the non-linear regime. At the smallest scale shown, k\ = 2.50 /iMpc -1 , the configuration 
dependence of Q is totally washed out. For the power spectrum, one-loop PT does remarkably well over the 
whole range of scales considered, remaining within less than 50% of the numerical simulation measurements 
and the fitting formulae up to scales where A(k) ps 100! 

Figure [)] shows the equilateral hierarchical amplitude Qeq for the n = —2, —1.5, —1 and CDM initial 
spectra as a function of scale. For the scale-free n = —2, —1.5 spectra, the one-loop predictions are those 
of equations ( |29] ) and ( |30] ) respectively. Note how well these simple formulae describe the behavior of Qeq 
from the tree-level value to the transition towards the strongly non- linear plateau. For the n = — 1 case, we 
do not show the one-loop correction, since as mentioned above, it does not obey self-similarity and therefore 
would not correspond to a single curve in the left bottom panel in Figure ^|. For the CDM case, we show 
the one-loop result for a$ = 0.33, which agrees very well with the numerical simulation measurements up to 
scales where A(A:) ss 10 (see left top panel in Fig. 0). The as = 0.64 output, however, corresponds mostly to 
scales already in the non- linear regime (see left top panel in Fig. ||); in this case equation (^0|) for equilateral 
configurations (not shown in Fig. ^) underestimates Qeq^), as expected from the discussion above. 
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An interesting question is whether the hierarchical amplitude Q actually becomes a constant independent 
of configuration in the highly non-linear regime. By comparing the results in Figure |] with the corresponding 
figures for the ki/k% = 2 configurations at the smallest scales, we see that the equilateral configurations 
seem to attain a slightly higher value for Q at large kR than the non-equilateral configurations. Given 
the uncertainties in our measurements, however, this trend is not statistically very significant, and more 
accurate measurements would be needed to assess in detail the validity of the hierarchical ansatz in the 
strongly non-linear regime. In fact, the simulations analyzed in this work do not probe very deeply into 
this regime. The validity of the hierarchical ansatz in the strongly non-linear regime has been considered in 
real space studies of the three and four-point correlation functions and counts in cells analyses (Bouchet & 
Hernquist 1992; Lahav et al. 1993; Colombi et al. 1994; Lucchin et al. 1994; Matsubara & Suto 1994; Suto 
& Matsubara 1994; Bonometto et al. 1995; Colombi et al. 1996; Ghigna et al. 1996.) After correction for 
finite volume effects, these results seem to show a small but significant departure from the scaling expected 
in the hierarchical ansatz. 



6. Conclusions and Discussion 

We have considered one-loop corrections to the bispectrum in Perturbation Theory (PT) and compared 
the results with numerical simulations for scale-free and CDM initial spectra with Gaussian initial conditions, 
focusing on the change of configuration dependence of the hierarchical amplitude Q as the transition to 
the non-linear regime is made. We found very good agreement between one-loop PT and our TV-body 
measurements for scale-free spectra with n = —2, —1.5, and for the CDM initial spectrum. For scale-free 
n > — 1 spectra, one-loop corrections diverge, and the simplest remedy of introducing a cutoff at small scales 
in the initial power spectrum (e.g., at the Nyquist frequency of the particles), breaks self-similarity, an effect 
which is not seen in the numerical simulations. On the other hand, for spectra with n > — 1, tree-level PT 
does well compared to numerical simulations even on scales comparable to the correlation length. 

For the power spectrum, we find excellent agreement between one-loop PT and the numerical simulations 
even on scales where the one-loop correction dominates over the tree-level contribution, i.e., where one would 
naively expect PT to break down. In fact, the simple expressions in equations (^) and ( |32| ) follow quite 
closely the fitting formulae for the non-linear power spectrum proposed by Jain, Mo & White (1995) and by 
Peacock & Dodds (1996) over a remarkable range of scales. For the hierarchical amplitude Q, we showed that 
one-loop corrections correctly describe the evolution of the configuration dependence observed in numerical 
simulations on weakly non- linear scales, for power spectra with sufficient relative large-scale power (n < — 1 
and CDM). At scales comparable to the correlation length, where one-loop contributions become of the 
same order as their tree- level counterparts, the numerical simulations show a progressive flattening of Q{0). 
This flattening starts at intermediate angles, as these configurations become increasingly probable due to 
"randomization" of density and velocity gradients, and propagates to collinear configurations (8 = 0,7r). 
One-loop PT does not reproduce this observed flattening very well, but it is nonetheless able to follow 
configurations close to collinear further into the non-linear regime. In the strong clustering regime, the N- 
body results show almost no dependence on configuration shape, Q ~ constant. This saturation value shows 
a clear dependence on the initial spectrum, consistent with previous numerical simulations in the literature 
(Efstathiou et al. 1988; Colombi et al. 1996), as parametrized by FMS, Q sat (n) m 3/(3 + n). Note that for 
the CDM case, a similar formula could be used, where n is taken as the effective spectral index at the scale 
of nonlincarity. In this case, we would find that, for the a% — 0.64 output time, n e s(k) « —1.8 at the scale 
where A(fc) = 1; this corresponds to Q S at(n e ft) ~ 2.5, in rough agreement with Figures ||| and ||. 
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In recent work related to our paper, Jing & Borner (1996) noted that the predictions of tree-level PT 
for the three-point function in real space did not agree with their numerical simulations of CDM models in 
the weakly non-linear regime. For this comparison, however, they considered scales close to the correlation 
length, for which one- loop corrections are expected to significantly alter the tree- level predictions. The 
results we obtain in Fourier space suggest that their measurements should agree very well with PT when 
one-loop corrections are included. 

Finally, in light of these results, we return to comment on the issue raised in the introduction, the shape 
of Q as a probe of bias and its robustness to non- linear effects. Although the non- gravitational effects that 
transform the non-linear density field into the observed distribution of luminous galaxies are undoubtedly 
complex, they may be relatively local; that is, suppose the probability of forming a luminous galaxy depends 
only on the underlying density field in its immediate vicinity. For simplicity, we also suppose this dependence 
is deterministic rather than stochastic. Under these simplifying assumptions, the relation between the galaxy 
density field S g (x) and the mass density field S(x) is of the form S g (x) = f(S(x)) = T] n b n 5 n , where b n are 
the bias parameters. For the reduced tree-level bispectrum, this local bias scheme implies 

Qf = ^Q (0) + | (33) 

(Fry & Gaztahaga 1993). Gaztaiiaga & Frieman (1994) have used the corresponding relation for the skewness 
S3 to infer b% ^ 1, 62 — from the APM catalog, but the results are degenerate due to the relative scale- 
independence of S3. Fry (1994b) has used the comparison between Q g inferred from the Lick catalog and 
the tree-level PT prediction Q^> to infer values for &i,&2. In particular, since Q g displays little of the 
configuration dependence expected of Q^°\ he finds a best fit value of 61 ~ 3, a large linear bias factor. On 
the other hand, in order to extract a statistically significant Q g (9) from the Lick catalog, an average over 
scales which include values of k comparable to the scale of nonlinearity was required. The question then 
becomes, to what degree do one-loop corrections to on these scales affect the determination of the b n 
from equation (J33])? 

To partially address this question, in Figure [lO] we show the expected correction to the hierarchical 
amplitude Q for k\jki = 2 configurations according to one-loop PT, for T = 0.25 CDM with erg = 0.64, 
at scales where Awl. Note that we do have TV-body results on these scales at this output time, but 
the statistical uncertainties are rather large due to the small number of independent modes available (see 
Appendix A.) The one- loop PT results in Figure [l^ agree within the errors with these measurements. For 
a realistic spectrum, Figure [To] illustrates how the configuration dependence of Q should change on scales 
relevant for observations. There is a minor although noticeable flattening of Q; if one were to instead attribute 
this to an effective bias according to equation (|3^), it would correspond to 1.25 < bi < 1.4, 0.5 < b 2 < 0.8, 
where low (high) values of bi are correlated with low (high) values of &2- This preliminary result shows 
that the estimate of the bias parameters from the measured bispectrum in the galaxy distribution using the 
tree-level result, equation (^), could be affected by non-linear evolution, unless scales much larger than the 
scale of non-linearity are considered. Further work is needed in order to quantify this issue better. 

In evaluating the prospects for measuring the bispectrum in galaxy surveys, observational considerations 
such as selection function, angular projection in two-dimensional surveys, redshift distortions, and aspects of 
survey design such as sky coverage, geometry, and sampling rate must be carefully considered. Although an 
exhaustive study of these kinds of questions is beyond the scope of this paper, we shall make some general 
comments. 

On large scales, for a scale- free power spectrum, P(k) ~ k n , the statistical uncertainty in Q scales with 
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configuration size as (k n g/kY 3+n ^ 2 per mode (FMS). To compare observations with perturbation theory, we 
need data on scales k < k n g and thus need a survey with many modes at the scale k n £, i.e., a survey that 
covers a large volume. In order to determine the optimal sampling strategy, it is necessary to know with 
precision what are the various sources of error, at least from the pure statistical point of view (e.g., Szapudi 
&; Colombi 1996). To reduce shot-noise uncertainties, it is desirable to construct surveys with a large number 
of galaxies, i.e., as complete as possible. In addition, minimizing edge effects requires compactness (i.e., the 
boundaries of the survey should have minimal surface), while finite- volume errors call for large sky coverage. 
The best sampling strategy results from balancing these various effects. Kaiser (1986) concluded that to 
measure the two-point function at large scales it is best to have sparse samples with large sky coverage and a 
sampling rate of order 1/10. In more recent work concerning the power spectrum, Heavens & Taylor (1997) 
reach similar conclusions. In the case where only a small part of the sky is covered, another issue arises: the 
choice of the catalog geometry. The conclusions of Kaiser (1996), who analyzed weak gravitational lensing 
statistics, favor a catalog made of many small patches spread over the sky. A similar study by Colombi, 
Szapudi, & Szalay (1997), based on counts-in-cells statistics (including higher order moments), reaches at 
least qualitatively the same conclusions; however, the latter conclude that the sampling rate should be 
increased as higher order statistics or smaller scales are considered. 

The discussion above can be illustrated by the situation in existing surveys, where there are on-going 
efforts to measure three-point statistics. In the QDOT and IRAS 1.2Jy surveys, the main limitation comes 
from discreteness effects, which dominate the signal even at large scales (Feldman et al. 1997). On the other 
hand, in the APM survey (Gaztanaga & Frieman 1997), the main challenge is to successfully deconvolve the 
angular projection, to extract the signal at quasilinear scales (Thomas & Fry 1997). The situation in the 
Las Campanas Redshift Survey is somewhat complicated, because of both geometry and selection function 
effects. In particular, the thin slices make the estimation of three-dimensional statistics rather uncertain 
(see, e.g., Heavens & Taylor 1997), even for the power spectrum case. For the bispectrum, the mixing of 
scales arising from the window function of a thin slice makes it at best difficult to separate out the quasilinear 
regime and compare with perturbation theory calculations. The prospects are much better for planned future 
surveys; in particular, the Sloan Digital Sky Survey and the Two Degree Field Survey should provide an 
accurate determination of the bispectrum over a wide range of scales in the weakly non-linear regime, with 
errors perhaps as small as a few percent in Q (Colombi, Szapudi, & Szalay 1997). 

Galaxy clustering derived from redshift surveys is distorted radially by peculiar velocities (Kaiser 1987). 
An important issue regarding the determination of bias using equation ([33] ) in redshift surveys is how redshift 
distortions are expected to modify the configuration dependence of Q. At large scales, tree-level PT predicts 
that redshift distortions increase the configuration dependence of for models with Q = 1 by 20%, whereas 
low f2 models show negligible redshift distortions (Hivon et al. 1995). This is as expected: peculiar velocities 
from coherent inflows, most important in the Q = 1 case, lead to more anisotropic structures in redshift space, 
thus increasing the configuration dependence of Q. On small scales, on the other hand, velocity dispersion 
makes Qeq less scale-dependent than in real space, as shown by modelling the velocity distribution function 
(Matsubara 1994) and numerical simulations (Lahav et al. 1993; Matsubara & Suto 1994; Suto & Matsubara 
1994; Bonometto et al. 1995; Ghigna et al. 1996). Although these works have concluded that higher-order 
statistics are therefore more hierarchical in redshift space than in real space, analysis of non-equilateral 
configurations leads in fact to the opposite conclusion: at small scales, the redshift space bispectrum shows 
increased configuration dependence due to anisotropics caused by velocity correlations along the line of sight, 
which enhance colinear configurations and suppress equilateral configurations relative to the real space case 
(Scoccimarro, Couchman & Frieman 1997). 
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An attractive feature of equation ( |33| ) as a tool to probe bias, is that the tree-level hierarchical amplitude 
Q(°) is very insensitive to the cosmological parameters Q and A, in contrast with determinations from 
large-scale flows which contain a degeneracy of the linear bias parameter bi with f2. It is interesting to 
check whether one-loop corrections to the bispectrum introduce a significant dependence on cosmological 
parameters. In Appendix B.3 we show that to a very good approximation, all the dependence of PT solutions 
on 17 and A can be described by the linear growth factor, to arbitrary order in PT. Therefore, for the same 
normalization of the linear power spectrum, or as , the hierarchical amplitude Q should be almost insensitive 
to Q and A even in the non-linear regime. We have checked this prediction against numerical simulations, and 
found that for the as = 0.64 output, the reduced bispectrum Q in an Q, = 0.5 (r = 0.25) open CDM model 
is virtually indistinguishable from the corresponding plots shown in Figure ||. The results in Appendix B.3 
also suggest that the same result regarding the f2 and A dependence should hold for higher order statistics, 
such as the S v parameters (p = 3, 4, . . .). Work is on progress on this issue (Jain & Colombi 1997). 

After this work was completed, we received a preprint by Matarrese, Verde, and Heavens (1997), which 
discusses error estimates on bispectrum measurements and a likelihood approach to extracting the bias. 
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A. Measuring the Bispectrum in Numerical Simulations 
A.l. The algorithm 

To measure the power spectrum and the bispectrum in the 256 3 scale-free simulations and in the CDM 
simulation, we wrote a FORTRAN program which, in brief, computes the density contrast on a grid by using 
"cloud-in-cell" (CIC) interpolation (see, e.g., Hockney & Eastwood 1988), fast Fourier transforms it, and 
then uses Monte-Carlo simulations for spatial averaging in fc-space. 

More specifically, given a triangle with sides k\, fe, and k^, the estimates of -P(fci), P(k2), P{k^) and 
B(ki, ^2, k$) are done as follows. The quantities we measure are actually smoothed over a bin of width Afc: 

i rk+Ak/2 
V{k) A-Afc/2 
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with 



and 



with 



P(q) = J %) S*(q) sin (9 d9 #, (A2) 

V(k) = 4irk 2 Ak[l + Afc 2 /(12fc 2 )], (A3) 

B(k 1 ,k 2 ,k 3 ) = 1 — — / B(qi,q 2 ,q 3 ) q 2 dqi q?,dq 2 qjdq 3 , (A4) 

1/ K2, K3J J qi e[fei-Afe,/2,fe,+Afci/2] 



„ 3 

B( qi ,q 2 ,q 3 ) = / Re^)^)^^)! TT(sin^^#0, (A5) 

iq 1 +q 2 +q 3 =o J £=1 

F(fei,fc2,fe 3 ) - / [S D }d 3 q 1 d 3 q 2 d 3 q 3 =8ir 2 k 1 k 2 k 3 Ak 1 Ak 2 Ak 3 . (A6) 

•/ o, 6 \ki -Aki In.ki +Aki Il\ 



In equation (A5), "Re" means "real part of". Most importantly, the function 6(q) is corrected for the effects 
of the smoothing in real space due to CIC interpolation, yielding 

^ = 8 sm(k x /2) sin( k y /2) sin(fc z /2) ^ 

where 8(q) is the actual Fourier transform of the density contrast computed on the grid. 

There is a constraint on the values of ki, k 2 , and k 3 so that the numbers qi € 2?i, where 

Vi = [h- Aki/2, h + Ah/2], (A8) 

form a triangle, that is 

k 3 > \ki — k 2 \ + 3Afc/2 and cyclic permutations, (A9) 

with 

Afc = (Aki + Ak 2 + Ak 3 )/3. (A10) 

As a result, if 9 represents the angle between vectors k\ and k 2 , the sampled values of 9 cannot be arbitrarily 
close to or close to tt: 

and 

r ffl , 3Afc|fc 1 -fc 2 |/2 + (3Afc) 2 /8 

[cos6»] max = 1 — . (A12) 



Now, let us imagine that we have chosen numbers hi, k 2 , and ^3 obeying the constraint of equation (At). 
Calculating integral (A4) is simply done by Monte-Carlo simulation, randomly choosing numbers qi in the 
intervals of equal probability T>i. We simultaneously compute integral ( A.5), i.e., we estimate the average of 
the quantity Ke[S(q 1 )S(q 2 )5(q 3 )] over all the possible positions of the solid body formed by the vectors qi, 
q 2 , and q 3 . The method used here consists of picking a random direction for q 1 and then choosing randomly 
the direction of q 2 in a circle of equal probability around q 1 such that the angle between q 2 and q 1 remains 
fixed (actually determined by the values of qi, q 2 , and q 3 ). Each iteration in our Monte-Carlo simulation 
thus consists of randomly choosing the numbers q%, q 2 , q 3 and the orientation of the solid formed by the 
vectors q l , q 2l and q 3 in space (three angles). 
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There is the problem that we have access to only discrete values of q on a three-dimensional grid Qij^k- 
This is solved by "random interpolation" . For a given value of q, we associate a probability to each of the 
eight nearest grid sites. These weights are computed the same way as for CIC interpolation. Each of these 
weights determines the probability of actually picking the corresponding value of q i j k , i.e., setting q = q i ■ k . 
To measure the bispectrum we apply this procedure to q 1 and q 2 an d se t q 3 = — <Zi — q 2 - 

Finally, one might wish to correct for discreteness effects, that is, to subtract off the contribution of the 
shot noise of the particles (see, e.g., Peebles 1980, eqs. [41.5] and [43.6]). We did so for all the measurements, 
except for the earliest output of the CDM simulation we analyzed, as explained in § 5. 



A. 2. Error Bars 

The method of estimating the errors in the figures depends on the case considered. For the scale-free 
PM simulations with 128 3 particles (n = —2, —1, 0, +1), since we have -/V rca = 4 independent realizations 
we can infer the errors from the dispersion of the measurements. If F is a quantity we measure, for which 
the estimator F is 

F = 



v rca _-, 



Nr, 

where Fi stands for the measurement of F in realization 1 < i < N rca , the estimator of the error reads 

[AF] 2 = — — — ^(-Fi — F) 2 . (A14) 

2—1 

For each of the two scale- free PM simulations with 256 3 particles, we had only a single realization to work 
with. In these cases, we extracted results at several output times and rescaled them under the assumption of 



self-similar evolution, forming the average of equation (A13). To estimate the errors, we treated the different 



output times as effectively different realizations and used the estimator (A14). However, the different output 
times are not actually statistically independent realizations. In these cases, the error bars on the figures are 
likely to be more a reflection of departures from self-similarity than of real underlying statistical uncertainties. 

For the CDM simulation, we have access to only one realization, and we cannot combine rescaled output 
times to artificially increase iV roa , because CDM is not scale-free. In this case, we use the error estimates 
of FMS and Feldman, Kaiser & Peacock (1994), which assume that the Fourier components are Gaussian- 
distributed: 

[AP(k)] 2 = J-[P tot (k)]\ (A15) 
V(k) 

[AB(k 1 ,k 2 ,k 3 )\ 2 = ^—^ [P t ot{ki)Ptot(k 2 )P tot (k 3 )} , (A16) 

2V(k ll k 2 ,k 3 ) 

where Ptot(fc) = P(k) + 1/A par is the total power. (We express wave numbers in units of the fundamental 
mode kf = 2tt/L.) Although the power spectrum and the bispectrum are statistically correlated, we use the 
standard error propagation formula to compute the error on the ratio Q(ki, k 2 , k 3 ). 



In equations ( Al5| ) and ( Al6| ), the quantities V(k) and V{ki, k 2 , k 3 ) represent, in units of kf, the 



"number of independent modes" for the power spectrum and the bispectrum. In principle, they should be 
equal to V(k) and V(ki, k 2 , k 3 ), but this does not make sense when the bins Afc and Aki are order of unity 
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or smaller (in units of kf), as is the case for our measurements (we chose Ak = Aki = 0.01). We are indeed 
working in a discrete Fourier space, in which the thinnest effective binning is Ak c g ~ 1. For Ak, Aki 1, 
we therefore take V(k) ~ 27rfc 2 and V r (fc 1 ,fc 2 ,fc 3 ) ~ (4/3)TT 2 k 1 k 2 k 3 , where we included symmetry factors of 



1/2 and 1/6 respectively (compare with Eqs. [A3] and |A6| ) . The reality constraint 6(h)* = S(—k) reduces 
the effective Fourier volume by one-half, whereas the three-fold symmetry of triangle configurations yields 
an additional factor of 1/3 (see FMS). 

In the 128 3 scale-free simulations, for which we have several independent realizations, we found that 



the errors one would estimate from equations ( A15 ) and (A16) are slightly smaller than the dispersion in 
equation (A14). This result suggests that equations (A15) and ( A16| ) underestimate the true errors. We 
do not know whether this is because of the Gaussian assumption quoted above or because the number of 
independent modes is overestimated. Rigorous calculation of the power spectrum and bispectrum error bars 
is a non-trivial issue that goes beyond the scope of this paper. 



In any case, neglecting the shot-noise contribution, Ptot(k) 
urations k\ = k 2 = k$, with the above assumptions one finds 



AQ 

Q 



l 



v/67rA(fc) 



P(k), and considering equilateral config- 



(A17) 



This means that the relative error on Q is expected to be larger on large length scales, which is unfortunate 
if one wants to probe the weakly nonlinear regime. Also, since one expects Q to decrease with spectral index 
n, the relative error on Q is expected to increase with n for the same degree of nonlinearity A(fc). 

Note finally that there is an error due to the finite number Niter of iterations used for the Monte-Carlo 
simulation discussed in the previous Section. A fair estimate of this error to first order is to use equations 
( |A15| ) and flAl6| ) with V(k) = N itcl or V{h,k 2 , k 3 ) = Niter- With our choice, Nit el = 10 7 , the corresponding 
uncertainty on the measurement of P(k) and B(ki, k 2 , k 3 ) is negligible compared to the other sources of 
error mentioned above. 



B. Eulerian Perturbation Theory 

B.l. The Equations of Motion 

Assuming the universe is dominated by pressureless dust (e.g., cold dark matter), in the single-stream 
approximation (prior to orbit crossing) one can adopt a fluid description of the cosmological TV-body prob- 
lem. In this limit, the relevant equations of motion in the Newtonian approximation to General Relativity 
correspond to conservation of mass and momentum and the Poisson equation (e.g., Peebles 1980, Scoccimarro 
& Frieman 1996). Assuming the initial velocity field is irrotational, the system can be described completely 
in terms of the density field and the velocity divergence, 9 = V • v. Defining the conformal time r = J dt/a 
and the conformal expansion rate H. = dlna/dr, the equations of motion in Fourier space become 



6{k,T) = - J d 3 fcx J d 3 k 2 S D {k - - k 2 )a(k,ki)e(k U T)S(k 2 ,r), (Bl) 



dd(k,r) 
dr 

+ H{r) 9(k, r) + ^nn 2 (r)~5(k, r) 



J d 3 h J d 3 k 2 5 D (k -ki- k 2 )[3{k, fej, k 2 )6{k x ,T)6{k 2 ,T), (B2) 
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where k is a comoving wave number, and 

a(k,ki) = ^-Jp-, 



/3(fe,fei,fe 2 ) 



fc 2 (fci -fc 2 ) 

2 J\i-y 



(B3) 



Equations ( JB l| ) and flB2| ) are valid in an arbitrary homogeneous and isotropic universe, which evolves ac 
cording to the Friedmann equations: 

dH(r) 



dr 



.^ 2 (r) + |a 2 (r), 
A 



(n-l)H 2 (r)=k--a 2 (T), 



(B4) 
(B5) 



where A is the cosmological constant, the spatial curvature constant k = —1,0, 1 for f2 tot < 1, fitot = 1) and 
fitot > 1, respectively, and fitot = fi + fiA, with £1a = Aa 2 /(3H 2 ). 



B.2. Perturbation Theory Solutions for 57 = 1 and A = 

For fi = 1, the perturbative growing mode solution for 8 is given by equation (^) and for the velocity 
divergence by 

oo 

9(k,T)=U(r)J2a n (r)e n (k). (B6) 
The nth order solution for S is given by equation ([}]), with a similar relation for the velocity field, 



d qi . . . id q„5i)(k - q 1 



l )G^{q l ,...,q n )8 1 {q 1 )...8 1 {q n ). 



(B7) 



The functions Fn an d Gn are constructed from the fundamental mode coupling functions a(k, fci) 
and (3(k,ki,k 2 ) by a recursive procedure (see Goroff et al. 1986; Jain & Bertschinger 1994), 



E 



Gm(<Zl, ' * * 5 Qm) 



G n (qt 



^ (2n + 3)(n- 1) 

+2/3(fe, k 1 ,k 2 )G n - m (q m+1 , ...,q 
G m (qi, ■ ■ ■ , <Z m ) 



(2n + l)a(k, ki)F n - m (q m+1 , q n ) 



) * * • ) In/ 



E 



3a(fc, fci)F n _ m (g m+1 , . . . , q n 



< t (2n + 3)(n- 1) 
+2n/3(k, k 1 ,k 2 )G n - m (q rn+1 , . . . , q„) 



(B8) 



(B9) 



(where fci = 
procedure: 



<Z,„, &2 = <Zm+l 



<7 n , fc = fci + fe2, and -Fx = C?x = 1), an d the symmetrization 



7r(ra)7> 



7T 



(BIO) 
(Bll) 



where the sum is taken over all the permutations 7r of the set {l,...,n}. Explicit expressions for the 
unsymmetrized kernels -F3 and F4 are given in Goroff et al. (1986). 
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B.3. £1 and A dependence of Perturbation Theory Kernels 

When Q 1 and/or A ^ 0, the PT solutions at each order become increasingly more complicated, due 
to the fact that growing modes at order n in PT do not scale as a n (r) as assumed in equations (||) and ( ]B(j| ) . 
Furthermore, the solutions at each order become non-separable functions of r and k (Bouchet et al. 1992, 
1995; Bernardeau 1994b; Catelan et al. 1995), so there appear to be no general recursion relations for the PT 
kernels in an arbitrary FRW cosmology. However, it is well known that the f2 and A dependence is extremely 
weak once the growth factors have been scaled out (Bouchet et al. 1992, 1995; Bernardeau 1994b). In this 
appendix we show that a simple approximation to the equations of motion leads to separable solutions to 
arbitrary order in PT and the same recursion relations as in the Einstein-de Sitter case. All the information 
on the dependence of the PT solutions on the cosmological parameters ft and A is then encoded in the linear 
growth factor, Di(t), which in turn corresponds to the normalization of the linear power spectrum, or ds. 



In linear PT, the solution to the equations of motion (Bl) and (JB2|) reads 



S(k,r) = D 1 (r)8 1 {k), (B12) 
9(k,T) = -H(T)f(n,A)D 1 (r)S 1 (k), (B13) 

where D\(t) is linear growing mode, which from the equations of motion must satisfy 
and /(Q, A) is defined as 

a ma n dr 

Explicit expressions or D\{t) and /(£1,A) are not needed for our purposes (see e.g. Peebles 1980). It 
is nevertheless important to note the following simple fits in the cosmologically interesting cases, namely 
f(fl, A) « fl 3 / 5 when A = (Peebles 1980), and /(fi, A) « ft 5 / 9 when fl + Q A = 1 (Bouchet et al. 1995). As 
mentioned before, we look for separable solutions of the form 



S(k,r) = J2D n (T)5 n (k), (B16) 

n=l 

oo 

6(k,r) = H(T)/(ft,A)^£J n (r)0 n (fc), (B17) 



n=l 



From the equations of motion (Bl) and (|B2j ) we get for the n order solutions 

n-l 



nf 



S n + E n e n = - J d 3 k 1 d 3 k 2 S D (k - fei - k 2 )a(k, fci) ^ D n ^ m E m 6 m (k 1 )S n - m (k 2 ), (B18) 



m— 1 



E 71 

—t 

nf 



3 n 



2/2 



3 n 



—r - l)E n 9„ + —^D n 5, 



2/2 



/ d 3 k 1 d 3 k 2 6 D (k - fei - k 2 )(3(k,k u k 2 ) E n - m E m m (ki)9 n -. m (k2), (B19) 

m=l 
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where the dot denotes a derivative with respect to r. By simple inspection, we see that if f(fl, A) = J7 1 / 2 , 
then the system of equations becomes indeed separable, with D n = E n = (_Di)™. In fact, the recursion 



relations then reduce to the standard Q — 1, A = case, shown in equations (B8) and flB9|). Then 
il// 2 = 1 leads to separability of the PT solutions to any order, generalizing what has been noted before 
in the case of second order PT by Martel & Freudling (1991). As mentioned above, the approximation 
/(0,A) rj fi 1 / 2 is actually very good in practice; for example, the exact solution for the A = case gives 
D2KD1) 2 = 1 + 3/17(Sl — 2 / 63 — 1), extremely insensitive to f2, even more than what the approximation 
/(0,A) — f2 3 / 5 w fi 1 / 2 would suggest, since for most of the time evolution fl and A are close to their 
Einstein-de Sitter values. 
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Fig. 1. — The left top panel shows the non-linear power spectrum in terms of A(fc) = Ank 3 P(k) as a function 
of scale for n = —2 scale-free initial conditions. Symbols denote measurements in numerical simulations, 
whereas lines show the linear, PD, JMW and one-loop perturbative results, as indicated. The other three 
panels show the hierarchical amplitude Q for triangle configurations with ki/k 2 — 2, as a function of the angle 
9 between fci and ^2, in numerical simulations and for tree-level and one-loop PT. The panels correspond 
to stages of non-linear evolution characterized by A(fci) = 0.71, 1.95, 3.92. 
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Fig. 2. — Same as Figure [l]for n=-1.5. 
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Fig. 3. — Same as Figure [j] and g for smaller (more non-linear) scales, showing the decrease in the configu- 
ration dependence of Q. The top panels correspond to the n = — 2 initial spectrum, the bottom panels show 
the n = —1.5 initial spectrum. 
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Fig. 4. — Same as Figure [l] for n — 
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Fig. 5. — Same as Figure [j]for n — (top panels), and n — 1 (bottom panels). 
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Fig. 6. — Same as Figure [l] for CDM initial power spectrum, with T = 0.25. These four panels correspond 
to a erg = 0.2057 output. 




Fig. 7. — Same as Figure ||for a later (more non-linear) output, as — 0.3291. 
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Fig. 8. — Same as Figure |]for the latest (most non-linear) output, erg = 0.64. 




Fig. 9. — The hierarchical amplitude Q for equilateral configurations as a function of scale for n = 
-2,-1.5, -1 and the CDM <r 8 = 0.33, 0.64 outputs. 



-34- 




Fig. 10. — The hierarchical amplitude Q for triangle configurations with ki/k 2 = 2 for the CDM cr 8 = 0.64 
output at scales where Awl. 



